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Abstract 

We apply the holonomic gradient method (HGM) introduced by [S] to the calcu- 
lation of orthant probabilities of multivariate normal distribution. The holonomic 
gradient method applied to orthant probabilities is found to be a variant of Plack- 
^ . ett's recurrence relation ([14]). However an implementation of the method yields 

i Q | recurrence relations more suitable for numerical computation than Plackett's re- 

currence relation. We derive some theoretical results on the holonomic system 
for the orthant probabilities. These results show that multivariate normal orthant 
| probabilities possess some remarkable properties from the viewpoint of holonomic 

CN ' systems. Finally we show that numerical performance of our method is comparable 

or superior compared to existing methods. 

Jsj ! 1 Introduction 

The holonomic gradient method (HGM) introduced by 
Nakayama et al. [9] is a new method of numerical calculation which utilizes algebraic 
properties of differential equations. This method has many applications in statistics. For 
example, an application to the evaluation of the exact distribution function of the largest 
root of a Wishart matrix is introduced in [5] and an application to the maximum likelihood 
estimation for the Fisher-Bingham distribution on the d- dimensional sphere is introduced 
in [7]. These applications greatly expand the scope of the field of algebraic statistics. 

In this paper, we utilize the holonomic gradient method for an accurate evaluation of 
the orthant probability 

$(£,/i) =P(X 1 >0,...,X d >0), (1) 

where the d- dimensional random vector X = (X\, . . . ,X^) € R d is normally distributed 
with mean /i and covariance matrix S, i.e., X ~ N(fi, S). Since evaluation of the orthant 
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probability has important applications in statistical practice, there are many studies about 
it. Genz introduced a method to calculate the orthant probability utilizing the quasi 
Monte-Carlo method in [4J. Miwa, Hayter and Kuriki proposed a recursive integration 
algorithm to evaluate the orthant probability in [S] . 

When the mean vector [i is equal to zero, the orthant probability can be interpreted 
as the area of a d — 1 dimensional spherical simplex (cf. [1]). In [16J, Schlaffi gave a 
classical differential recurrence formula for $(E,0). Plackett generalized Schlaffi's result 
and gave a recurrence formula for $(E,/z) in [14] . He evaluated orthant probabilities 
by a recursive integration utilizing his formula. Gassmann implemented the Plackett 's 
method in the case of higher dimensions in [3] . Their method is a recursive integration 
based on a differential recurrence formula, whereas the holonomic gradient method utilizes 
differential equations. Plackett's recurrence formula is not suitable for the holonomic 
gradient method and in this paper we give a new recurrence formula, which is more 
natural from the viewpoint of holonomic gradient method. 

Let y = and x = — jS -1 . We denote the i-th element of y by yi and the 

element of x by Xij. By this transformation of the parameters, the orthant probability in 
([I]) can be written as 



det(x) 1/2 exp ^-^x r yJ g(x, y), 



where 

g(x, y) = I ■■■ I exp I \ ] Xijtitj + y ]y{ti 1 dt (dt = dt l ... dt d ). (2) 



/oo p oo / d d \ 

. . . / exp I ^ x iMj + ^2 Viti I 
— 1 i — 1 / 



In order to evaluate the orthant probability, it is enough to evaluate g(x,y). 

To apply the holonomic gradient method, we need an explicit form of a Pfaffian system 
(0, Section 2]) associated with g(x,y). The Pfaffian system can be obtained from a 
holonomic system for g(x,y) (see [9]). For a given d, we can obtain a holonomic system 
for g(x,y) by applying the algorithms introduced in [T2] with computer algebra systems. 
However, for general d, these algorithms can not be applied and we need theoretical 
considerations. In this paper, we give a holonomic system for the function g(x,y) and 
construct a Pfaffian system from the holonomic system. 

Note that the integral g{x,y) satisfies an incomplete A-hypergeometric system QlOj) 
when the integration domain is not the orthant but a polytope. 

The organization of this paper is as follows. In Section [2J we describe a holonomic 
system for an integral 

g(x,y) = I f(t)exp ( ^ XijUtj + S Ty i ti ) dt, (3) 



. / d d \ 

' f(t) exp S2 x ijtitj + ^ yiU dt, 
Rd W=i i=i / 



where f(t) is a holonomic function or a holonomic distribution. In Section |3l we construct 
a holonomic system for the function in ([2]) by utilizing the result in Section [2j Then we 
construct a Pfaffian system from the holonomic system. Finally in Section H] we describe 
numerical experiments of the holonomic gradient method. 
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2 Holonomic system associated with the expectation 
under multivariate normal distributions 



In this section, we consider the holonomic ideal which annihilates the integral ([3]), which 
is the expectation under multivariate normal distributions (for the definition of the holo- 
nomic ideal, see |15j). We develop a general theory, where f(t) in ([3]) is a smooth function 
or a distribution in the sense of Schwartz (|17J). This theory is a generalization of the re- 
sult introduced in [5J. In Section[3]we will specialize fit) to be the indicator function of the 
positive orthant. Note that the results of this section can be applied to various problems 
of the multivariate normal distribution theory, other than the orthant probability 

We denote the ring of differential operators in x with polynomial coefficient by D = 
C(x\, . . . , x n , d Xl , . . . , d Xn ). The operator d/dxi is denoted by d Xi . We frequently use the 
following rings: 

D xyt := C(x {j , y k , t k , d Xij ,d Vk ,d tk : 1 < i < j < d, 1 < k < d), 
D xy := C(xij, y k , d Xij ,d yk : 1 < i < j < d, 1 < k < d), 
D x := C{xij, d Xij :l<i<j<d), 

D t :=C(t h d u :l<i<d). (4) 
We use the following notation. 

^ d d 

y = E"V> x = --S _1 , h(x, ^, t) = Xjjtjtj + yjtj. (5) 

i,j=l i=l 

2.1 The case of a smooth function 

We first consider the case when f(t) is a smooth function. In this case, the integral ([3]) 
converges if the matrix — x is positive definite and f(t) is of exponential growth. 

To state the main theorem of this section, we show a general lemma. The notation 
/, g, x, y in the following lemma is generic and not related to ([5]). 

Lemma 1. Let D x (resp. D xy ) be the ring of differential operators with polynomial coeffi- 
cient C(xi, d Xi : 1 < i < n) (resp. C(x{, yj, d Xi , d yj : 1 < i < n, 1 < j < m) ). Suppose that 
a holonomic ideal I annihilates a function f(x, y) on TV 1 x R m and the function f(x, y) is 
rapidly decreasing with respect to the variable y for any x in an open set O C R n . Then 
the integration ideal of I with respect to the variable y annihilates 

9i x ) / f(x,y)dy (6) 

JR, m 

defined on the open set O. 

Proof. Since the function f(x,y) is rapidly decreasing, the integral of ([6J converges. Let 
P be in the integration ideal J = (i + Y^T=i ®VjD x> y )^D X , then we have P»g = J P»fdy 
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by the Lebesgue convergence theorem. Since the differential operator P can be written 

as 

m 

p = p Q + ^d Vi P i eD x (P Q ei, PieD xy ), 

we have J P • fdy = ^ J Vi Pi • fdy. Since Pi • f is rapidly decreasing, we have J d Vi P{ • 
fdy = 0. ^ ' □ 

We now go back to g(x,y) in ([3]) and use the notation in (j3J) and (jSJ). Consider a 
C-algebra morphism from D t to D^j defined by 



<p : A -> D 



Here, we assume x^ = x ji} Xij = X]V Since [tp(d ti ), <p{d tj )\ := ip(d ti )ip(t j ) - tp(tj)(p(d ti ) = 
5ij, where 5ij is Kronecker's delta, ip is well-defined as a morphism of C-algebra. 
Now we state the main theorem in this section. 

Theorem 2. Suppose a function f on R d is smooth and of exponential growth. If dif- 
ferential operators Pi, . . . , P s G D t annihilate f , then the following differential operators 
annihilate the integral g(x,y) in (J3J). 

<p(Pk) (l<k< s), (7) 
d Xij -2d yi d y . (l<i<j<d), (8) 

d*u-% (l<<<d). (9) 

Moreover, the differential operators (|7j), (jHJ), © generate a holonomic ideal in D xy if the 
differential operators P\, . . . ,P S generate a holonomic ideal in D t . 

Proof. For a differential operator 
and pi, q { G (i = 1, . . . , d), we put 

m ft) = E • • -^i 1 • • • ( 10 ) 

By the assumption, the differential operators 

P e {£=l,...,s), d Xi . (l<i<j<d), d yi {l<i<d) 

annihilate / and generate a holonomic ideal in D xyt . By the lemma introduced in [131 
Section 3.3], the differential operators 

Oh Oh Oh 

Pe(U,O k -—)(£=!,. ..,s), d Xij - — (l<i<j<d), d m - — {l<i<d) 



annihilate exp(h)f and generate holonomic ideal I in D xyt . As we will show in Lemma 

H] below, the differential operators (J7J), ©, © generate the integration ideal J of I with 
respect to the variables £j (i = 1, . . . , d). 

Since the function exp (h) f is rapidly decreasing, the integration ideal J annihilates 
the integral g(x,y) by Lemma dJ Hence, the differential operators fl7D,flSJ), 
© annihilate g(x,y). 

Since the integration ideal of a holonomic ideal is also holonomic (see, e.g., [2j Chap 

I] ), the ideal J is holonomic when the ideal / is holonomic. □ 

Now, we show that the integration ideal J is generated by the differential operators 
© , © , © in the following two lemmas. We denote P = Q when P — Qe J2t=i D xy t{d yi — 
/,)• 

Lemma 3. Let p t = d u - y» - ^Yfk=i x ikU and & = d u - y% - ^Yl=i x ikd Vk - For a 
differential operator P G D t and a multi index a G Nq = {0, 1, 2 . . . } d , the following 
equivalence relations hold. 

PfaPi) = P(d yi ;qi), (11) 
t a P(d y -q t ) = d«P{d y - qi ). (12) 



Here, we use the notation in flTO]) . 

Proof. By the straightforward calculation, we have the following relations for 1 < i, j < d. 

\Pi,Qj]=0, (13) 

\Pi,Pj] = [Qi,Qj]=Q, ( 14 ) 
[q u t j } = [q i ,d yj ]=5 ij . (15) 

In order to prove ( TTTT) . it is sufficient to show 

t? ■ ■ ■ t a M x = a* • • • d£q? 1 • • • q^ (a, (3 G Nq). (16) 

By the induction on we have a relation 

t^^dy-q* (l<i,j<d, AG N ). (17) 

When (3i = 0, the relation ffTTl) holds clearly. Suppose that the relation (TTTT) holds for 
Then we have 



■ Ji — i 



(qid yj - <%) = 9 w gigf* = ^.gf +1 . 



By induction, the relation ( ITTj) holds for any 
By the relations (fl3l) and ( TT71) . we have 



= (AeN ). 
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Now we prove the relation f lT6|) by induction on the multi index a G Nq. When a = 0, 
the relation ( TIB]) holds because of f|T8|) . Suppose the relation f fl6|) holds for a, then we 
have 

+ .+ ai . . .i ad r,^ ■ ■ ■ rP d — i f) ai ■ ■ ■ F) ad n^ 1 ■ ■ ■ n^ d 
L i L l l d Pi Pd — L i u yi °y d Vl Id 

= d m ■ ■ - d 2 (II ^ ® and (USD) 

= W • • • • • • $ dSl) and (PD). 

By induction, the relation (TlTj) holds for any a. 

Finally, the relation ([12} holds bv (III J) . ([15]) and (ITTjl . □ 

Lemma 4. T/ie integration ideal J is generated by the differential operators 
©,©,©■ 

Proof. Note that the ideal / is generated by differential operators 

d 

Pt{U]d u -yi -2^2x ik ti) (£=l,...,s), 

k=i 

0, i; --ltJj (l<i<j<d), 

d Xu -tl dn-U (l<i<d). 

By (|TT]) . the ideal / is generated by 



d 



P i (d yi ;d ti -y i -2j2xikd yk ) {l=l,...,s), (19) 

k=i 

d Xij -2d yi d yj (l<i<j<d), (20) 

d Xii -dl (1 <*<£*), (21) 

d yi -U (l<i<d). (22) 

We denote by J the left ideal generated by (0), (EJ), ©. Clearly we have J C J. If P 
is a differential operator in J, then P can be written as 

d d 

p = Q + J2 9 ti Ri + J2 S * ^ - *<) e St E D xyt ) (23) 

2=1 1=1 

by the definition of integration ideal. The differential operator Q is written as a linear 
combination of the differential operators ([T9~[) - f[2~T[) with D xyt coefficients. By the second 
term of the right-hand side of ([231) . we can assume without loss of generality that the 
variables d ti do not appear in these coefficients. By the relation f[T2"j) in Lemma [21 we can 



6 



assume that the coefficient in Q is an element in D xy . Since any differential operator in 
( TT9|) ( ]2Tj) is an element of J + ^2 i=1 d ti ■ D xyt , we can assume Q e J. 
Consider the equation 

d d 

P-Q- y £ i S i {a yi -t i ) = ^2d tl R t . 

i=l i=l 

We can assume that the variables d tl , . . . , d td do not appear in Si, . . . , Sj. For example, if 
tid tl is a term of S 2 then we can replace S 2 and Ri to S 2 — tid tl — 1 and Ri — ti(d y , 2 — t 2 ) 
respectively. The term tid tl in S 2 is removed. In the same way, we can remove all terms 
which include d+ . 

Expanding both sides and comparing the coefficients of dt t , we have 

P-Q = Y / Si(d yi -t i ). 

The right-hand side of this equation is an element of the left ideal J' := Yli=i D xy r (d Vi —U). 
Let the weight of tj be 1 and that of other variables 0, and consider a term order -< with 
this weight. The set {i; — d m |1 < i < d} is a Grobner basis of J' with the order, so that 
the initial term of P — Q has to divide some £«. Since P — Q is in D xy , we have P — Q = 0. 
Thus P E J. D 



2.2 The case of a non-smooth function 

Next we consider the case when f(t) is not smooth. In this case we consider f(t) as a 
distribution in the sense of Schwartz ([T7]). 
Let Q be a domain defined by 

{(x,y)\— x is positive definite}. 

For a tempered distribution / on R d , we can define a function on Q as 

g(x,y) = {f,exp(h(t,y,x))). (24) 

Since exp (h(t, y, x)) is rapidly decreasing with respect to the variable t when — x is positive 
definite, the right-hand side of ( )24l) is finite. 

A holonomic system for ( 1241) is given as follows. 

Theorem 5. // differential operators Pi, . . . ,P S G D t annihilate a tempered distribution 
f on H d , then the differential operators ([7]), OH]), © annihilate the function g(x,y) in f f24|) . 
Moreover, the differential operators (Ej), (IH]), (EJ) generate a holonomic ideal in D xy if the 
differential operators Pi,...,P s generate a holonomic ideal in D t . 

Proof. We only need to prove that the differential operators (J7|), (JSJ), (JSJ) annihilate g(x, y). 
Let P e = Yl c aP t a d^ and P/ = £ c a ^d?. We have 

d 

p e(d yi ; -yi-2^2 x ikd yk )(f, exp (h(t, y,x))) 

k=i 
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d 



(/, Pi{d Vi ;-yi-2^2 x ik d yk ) exp (h(t, y,x))) 



k=i 



(f, Pi{9 yi ] ~d u ) exp (h{t, y, x))) 
(f, P;(-d ti ;d yi ) exp (fc(t, y, x))) 
(/^/(-d^exp^t^x))) 
(Pt(ti;d t .)f,ex.p(h(t,y,x))) 
0. 



□ 



3 Holonomic system associated with the orthant prob- 
ability 

In this section we specialize fit) of the last section to the indicator function of the positive 
orthant and we will construct a Pfaffian system associated with the integral ([2]) for the 
orthant probability. 

3.1 Generators of the holonomic ideal 

At first, we obtain generators of a holonomic ideal which annihilates ([2]) by Theorem [51 
Let E be the positive orthant in R d defined by 



and Is be the indicator function of E. A holonomic ideal which annihilates 1e is given 
as follows. 

Lemma 6. The indicator function 1e is annihilated by the following differential operators 
as a distribution. 



The differential operators ( |25l) generate a holonomic ideal J in the ring D t . 

Proof. At first, we show that the differential operators Ud^ annihilates the function 1^. 
It suffices to prove for i = 1. Let cp(t) be a rapidly decreasing function on R d . Then, we 
have 



for any tj G R, 2 < j < d. Integrating both sides with respect to the variables t 2} . . . , t d , 
we have {tid tl lE,¥>) = 0. Therefore, the distribution tid tl lE is equal to 0. 

Next, we show that the left ideal J is holonomic. By the Buchberger's criterion, we 
can show that the set of the differential operators in (T25]) is a Grobner basis of J with the 



{t=(t u ...,t d ) eR% >o,(i 



l,...,d)} 



tid tl , . . . , tdd td 



(25) 
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weight w = (0, 1). The characteristic variety (see, e.g. [TTJ, [15]) of J is ch(J) = {(£,£) £ 
C 2c( |tj^ = 0, 1 < i < d}. The variety ch(J) can be decomposed as follows, 

U {(t,Z)\t i = o,s j = o,ieJ,j<£J}. 

Jc{i,...,4 

Since the Krull dimension of the each component is d, we have dim(ch(J)) = d and the 
ideal J is holonomic. □ 

The function g(x, y) in (|2J) can be written as g(x, y) = (le(t), exp y, £))). This is a 
case of Theorem [5] in which the distribution / is 1^. A holonomic ideal which annihilates 
1e is given in Lemma [61 Hence we have the following theorem for g(x,y) in (J2J). 

Theorem 7. Differential operators 
d 

2 x ikd Vl d yk + y<9 tt + 1 (i = 1, . . . , d, x {j = x Si ), (26) 

k=l 

d x ..-2d yi d V] (l<i<j<d), (27) 
d Xii -dl {l<i<d) (28) 

annihilate the function g{x,y) in (T5]) ; and generate a holonomic ideal I in the ring D xy . 
3.2 Differential recurrence formula 

Next, we give a Pfaffian system associated with (j2J). In order to write the Pfaffian system, 
we define new notation. For J C [d] — {1, . . . , d}, we put 

(x, y, x) — Xijtitj + ?/fc£fc, 

ieJ jeJ fcgj 

poo poo 

9j(x,y) = ... j exp(hj(x,y,t))dtj, 
Jo Jo 

where dtj = Yljej dtj. When the set J is empty, we set g$ — 1. For example, the functions 
are written as follows when d = 2. 







poo 


9{1,2}(X; 


y) 


= / exp 






Jo 






= g(x,y), 






POO 


9{1}{X- 


y) 


= / exp 






Jo 






POO 


9{2}(x, 


y) 


= / exp 






Jo 






= 1. 
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(30) 



Let J be a subset of [d] and put 

Qj = -^i + 2|>Aj (j = l,...,d), (29) 

Qj,j = - ( y 3 ■• + 2 x Jkd yk {j e J). 
\ kej J 

Lemma 8. The following equation holds. 

QjQj = Qj,jgj = 9j\{j} U e J). (31) 

Proof. Since gj is constant with respect to ye (£ ^ J), we have Qjgj = Qjjgj- We assume 
that J = [d] without loss of generality. Applying Qj to the integrand of gw, we have 

Vi + 2 ) exp W 37 ' ^' *)) 

fc=i / 

yj + 2 ^ rE jfe tfc ) exp (/i(x, £)) 
v fc=i / 
= -<9 % exp (h(x,y,t)) . 

Integrating both sides of the equation from to oo with respect to tj, we have 

/■oo 

/ exp(h(x, y, t))dtj = exp (h [dl \ {j} (x, y, t)) . 
./o 

Integrating both sides of the equation with respect to the remaining variables, we have 
the equation ( 131]) . □ 

Remark. By Lemma [BJ we have QiQjg = 9[d\\{i,j} f° r i 7^ j When the vector y is 
equal to zero, this equation can be written as 

Xi k Xjed Xke g = g[d\\{i,j}- 

k=l l<k^£<d J 

By the transformation of parameters with £ = (cr^) = — we have 

9 1 1 , 1/2 N I 1—1/2 

— T ' g) = \a\ 1 g [d ]\{i,j}. 

It can be easily checked that this equation corresponds to the classical Schlaffi's formula. 
For J = {ji, . . . ,j s } C [d], (ji < j 2 < • • • < j s ), we denote 
X J = ( x j k jJi<k,e<s, yj = (y jl ,...,y js ) T , 

= ~2 X f = VJ = ^jyj= (nj v ...,nj s ). 

The following differential recurrence formula holds. 
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Theorem 9. For any J C [d]. 

d = (^ + E 3eJ 4Mw J (32) 

{2d yi d yj gj {i,j}cJ,i<j 
d^gj {i} C J, i = j (33) 
else. 

Proof. For i J, we have c^gj = since the gj is constant with respect to i/i. For « 6 J 
Lemma [S] implies 

/^#J + a ij9j\{j} = \49J ~ Y a i ( y J i + 2 S X /A J # J 

= [ri - J2 "iiVs ~ 2 Y1 a i x id yk I gj- 

We have a ijVj = A*/ by the relation /i J = £ J ?/ J and we also have 

-2 (T^ k d Vh = Y 5ikd y* = 9 w> 

k,jeJ fceJ 

since — TL 3 x J is the identity matrix of size \J\. Hence, the right-hand side of (132]) equals 
d Vi gj- 

For {«, j} ^ J, d Xij gj = since the gj is constant with respect to x^. For {i, j} C J 
9^ exp y, t)) = (2 - 5 ij )d y% dy ] exp y, t)) . 

Integrating both sides we have the relation (13"31 . □ 

By theorem [9], we have differential operators d Xij — Ay, <9 yi — Aj which annihilate the 
vector value function G(x,y) = {gj{%,y))jc[d\- Here, Ay and Aj are 2 d x 2 d matrices 
with rational function entries. In the next subsection, we prove that the system of these 
differential operators is a Pfaffian system in the meaning of [TJ Section 2]. Note that the 
Pfaffian system have no singular point on 

{(x,y)\— x is positive definite }. (34) 



3.3 Pfaffian system and the holonomic rank 

In this section, we give the holonomic rank of the ideal I generated by (|26l) -(f28l). and 
show that the differential recurrence formula ( 132]) and (133]) in the subsection 13.21 give a 
Pfaffian system associated with (J2]). 
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In this subsection, we denote the ring of differential operators in the variables x, y by 
R. The holonomic rank of I is the dimension of R/RI as a vector space over the field of 
rational functions C(x, y) (see, e.g., [9]). 

At first, we give a lower bound of the holonomic rank. Since the holonomic rank 
equals to the dimension of holomorphic solutions of I • / = at generic points, we can 
obtain a lower bound of the holonomic rank by constructing linearly independent functions 
annihilated by /. 

Lemma 10. The holonomic rank of the ideal I in Theorem^ is not less than 2 d , i.e., 
rank I >2 d . 

Proof. For a vector e — (ex, . . . , e^) G {±l} d , let E £ be an orthant 

{t = (tx, ...,t d )e R d M, > (i = l, . . . , d)} . 

It is enough to show that the following 2 d functions are linearly independent and annihi- 
lated by J; 

9e(x,y)= exp(h(x,y,t))dt. 
Je £ 

By an analogous way as in the proof of Lemma El we can show that the indicator 
function of E £ is annihilated by the differential operators in (125]) for any e G {±l} d . 
Analogously to the proof of Theorem we can show that the differential operators fl26]) - 
f[2"gj) annihilate g e (x,y). 

Let c £ be a real number for e G {±l} d , and suppose ^ e c £ g £ = 0. Multiplying both 
sides of the equation by (27r)~ d//2 (det S)n _1//2 exp(— i/ZE" 1 ^), we have 

£ c e P(£ e |/i,£) = 0. 

e£{±l} d 

Here, P{E E \fi, E) is the probability of the event E e under the multivariate normal dis- 
tribution N(fjL,Ti). Substituting [i = te,(e G {±l} d ) and taking a limit t — > +oo, we 
have 

c £ = 

Hence, the functions g E (x,y) are linearly independent. □ 

In order to obtain an upper bound of the holonomic rank of /, we construct bases of 
R/RI as a linear space over C(x,y). The bases correspond to the functions gj. 
For J C [d], we put a differential operator 

pj= n Si'* ( 3s ) 

j'e[d]\J 

where is the differential operator in (j29|) . Note that the differential operators Qj 
commute with each other. By Lemma [H we have 

Pj9 = 9 J- (36) 
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Equation ( 136]) means that the differential operator Pj corresponds to the function g,j. For 
example, When d — 2 and J = 0, the equation fl36|) is written as follows. 

(y 1 + 2x u d yi + 2x 12 d y2 ) (y 2 + 2x 2 id yi + 2x 22 d y2 ) g { i >2 } = 1. 

Since the differential operator Qj commutes with d Xij — 2d yi d Vj (1 < i < j < d) and 
d Xii — dy. (1 < i < d), we have the following lemma. 

Lemma 11. The following formulas hold in R/RI. 

d Xij Pj = 2d m yj Pj (1 < i < j < d, J C [d]) , (37) 
d Xii Pj = d 2 y Pj (l<i<d,J C[d]). (38) 

Proof. For 1 < i < j < d and 1 < k < d, we have 

d^Qfc = ~d Xij ^y k + 2^2x ki d y ^j = Q k d Xij - 25 ik d Vj - 2S jk d yi , 

2d y% d yj Q k = -2d m d yj ^y k + 2 s ^x ke .d y }j = Q k 2d yi d y] - 28 ik d y . - 25 jk d m . 



For 1 < i < d and 1 < k < d, we have 
d Xii Q k = ~d Xii y k + 



2 / J Xkjdyt I — Qkd Xll — 2S ik d yi , 
i=i J 



d lQk = ~ 9 l V* + 2 ^2 X ^ d Ve = Qkd yi - 25^. 



^=1 



Therefore, the differential operator Qj commutes with d Xi . — 2d yi d yj (1 < i < j < d) and 
d Xii — dy. (1 < % < d). Then we have 

(d Xij - 2d yi d yj )Pj = Pj{d Xij - 2<)„ i) !h ) = 

in R/RL Similarly we have (d Xii — dy.)Pj = 0. □ 

The following lemma corresponds to Lemma [H 
Lemma 12. With the same notations as in Lemma El 

Q 3 Pj = Q jtJ Pj = P A{J} (j G J, J C [d]). (39) 

holds in R/RL 
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Proof. By definition of Pj, it is clear that QjPj = Pj\{j}- Since d ye Qe G /, dgPj is in I if 
£ J. Hence we have 



QjPj = - f yj + 2 x ikd Vk J JJ Q£ = Q 
\ k=i J ie\d\\j 



hJ r J 

[d]\j 

in R/RL Note that Qj and commute if j ^ k. □ 
Theorem 13. The following relations hold in R/RI for any J C [d\. 

d yi Pj = |^ + ^^^ [ eJ (40) 

(2d yi d yj Pj {i,j}cJ,i<j 

d ^ p J = Wn P J {i}cJ,i = j (41) 
[0 e/se. 

Proof. For z G J, we can show the equation ([40]) by Lemma [12] and an analogous calcu- 
lation as in the proof of fl32|) . For i ^ J, the equation ( 140]) is shown as in the proof of 
Lemma IT21 By Lemma [TT| we have the equation fT4"T|) . □ 

Corollary 14. T/ie se£ of the differential operators Pj (J C [<i]) in ( 1351) spans the quotient 
space R/RI as a vector space over C(x,y). 

This corollary together with Lemma [TU] establishes the following theorem. 

Theorem 15. 



rank 1 = 2 



d 



4 Numerical experiments 

In this section we present numerical experiments of our holonomic gradient method for 
orthant probabilities. Our experiments show that the holonomic gradient method is very 
accurate and fast compared to existing methods. 

By theorem [9] we have an explicit formula for ^G(x(t),y(t)) when x(t) and y(t) are 
smooth functions. In order to evaluate G(x,y) at (xi,yi), we put 

x(t) = (1 -t)x + txi, y(t) =ty x < t < 1. (42) 

Here, Xq is the diagonal matrix whose (i,i)-entry equals to that of X\. The initial value 
can be written as 

^(0)^(0)) = IJ (-T^iC )) ^ ■ ( 43 ) 
jeJ 

Note that j^G(x(t),y(t)) does not have singular points on [0, 1] since the Pfaffian system 
for G(x,y) does not have singular point on (]3"4"]1. 
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Table 1: Errors 



IN O. 


dim 




1 

L 


o 




2 


3 


5.473549e-08 


3 


4 


3.373671e-08 


4 


5 


2.265284e-09 


5 


6 


1.120033e-08 


6 


7 


7.330036e-09 


7 


8 


8.705609e-09 


8 


9 


2.288549e-09 


9 


10 


5.024879e-10 



The accuracy of the holonomic gradient method can be checked by looking at the 
summation Y^ £ e{±l} d P(-X" G £J e ) = 1. Table [TJ shows errors |1 — J2 £ e{±i} d e -^e)l f° r 
sample data. 

For a correlation matrix R = {pij} with pjj = p, z ^ j, the orthant probability 
can be written as a one-dimensional integral. Table [2] shows the values of the orthant 
probability calculated by the one- dimensional integral and HGM for certain values of p. 
The differences of the result of HGM and the one-dimensional integral are less than 10 -06 . 



Table 2: Comparing of one-dimensional integral 



rho 


Dunnett 


HGM 


error 


0.0 


0.00097656249807343551 


0.000976562500000 


1.926564e-12 


0.1 


0.0065864743711607976 


0.006586475124405 


7.532442e-10 


0.25 


0.026603192349344017 


0.026603192572268 


2.229240e-10 


0.5 


0.090909089922689604 


0.090909086078297 


3.844393e-09 



Table |3] shows averages of computational times of holonomic gradient method and 
Miwa's method for 100 sample data. The mean vectors of the sample data are all zero, 
and the covariance matrices are randomly generated. Table [3] shows that our holonomic 
gradient method evaluates orthant probabilities faster than Miwa's method, as the di- 
mension becomes larger. 

However, it should be noted that the holonomic gradient method can be very slow 
when the mean vector is far away from zero. When the mean vector is far away from zero 
and some eigenvalue of the covariance matrix is very small, the value of the integral ([2]) 
becomes very large since the parameter y becomes large. In such a case, the Runge-Kutta 
method takes a long time. 
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Table 3: Averages of computational times 



dim 


Miwa 


HGM 


dim 


Miwa 


HGM 


5 


0.002 


0.016 


9 


6.078 


1.050 


6 


0.011 


0.056 


10 


60.171 


2.371 


7 


0.080 


0.154 


11 


671.370 


5.411 


8 


0.664 


0.390 


12 




13.48 
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